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I.  It troduction 


Various  optimum  and  sub-optimum  filters  have  been  added  to  the 
basic  generlized  cross-correlator  in  order  to  enhance  its  detection  and 
estimation  capabilities  in  the  present  of  noise.  An  experimental  compar- 
ision  of  the  optimum  and  sub-optimum  filters  was  discussed  in  Hassab  and 
Boucher £1  ]  .For  the  random  signals,  the  optimum  window  which  they 

derived  has  the  results  of  time-delay  estimation  considerably. 

Basically,  the  solution  consists  of  cross-correlatin..;  the 
sensor  outputs.  The  time  argument  that  corresponds  to  the  maximum  peak 
in  the  output  is  the  time  delay(Fig.  l).  In  this  report,  time  delay 
estimation  under  several  assumptions  about  signals  is  obtained  by  averaging 
the  periodgrams  of  the  data  segments^ When  the  signal  is  unknown,  the 
maximum  entropy  power  spectrum  estimation  of  Burg's  algorithm  is  used  to 
estimate  the  signal  power  spectrum  from  the  measurement  of  separated 
sensors.  The  Ak-uke  FPK(  Final  power  error)criterion  is  also  concerned  in 
deciding  the  order  of  AH  model  and  signal  &  noise  power  spectrum  estimation 
of  Burg's  algorithm. 

TI.  Theory 


Let  the  sei.sor  output  be  the  form 
Z^(t)=y(t)+n ((t) 

Zp(t)=ay(t-T)+n2(t)  (l) 

where  the  signal  ,y(t)  and  noises  n^(t),n2(t)are  uncorrelated  and  jointly 
stationary  zero  mean  Qaussian  random  processes  over  an  observation'  interval 
>>1*.  Here  Tie  the  difference  in  arrival  time  or  the  time  delay  between  the 
two  sensor  elements. 

In  the  correlator,  the  multiplier  output  is  given  by(2} 

Z1*(w)Z;?(w)  =  a  4>y(w)  e“jW  +  4n1n?( “4^ w)  (£> 


where 


where 

4>y(w)*|Y(w)|2  ,  ♦n1n?(w),sN1  (w)  N?(w) 

(w)=Y*(w)N  (w)+aY(w)e  ^WrN  *(w) 

C  u  ! 


The  lower  case  <f>  refers  to  a  single  realization  of  the . function  and  the  upper 
case  indicates  its  expected  value  over  the  ensemble,  i.e.<  0y(w)>rr<  |Y<w)|‘>=  ^(w 
The  window  W^(w)  derived  by  Hassab  and  Boucher  £  p  "]  is  based  on  the  optimal 
detection  criterion  of  maximizing  the  ratio  of  mean  correlator  output  due  to 
signal  and  noise  present.  The  window  function  can  be  expressed  in  the  simple 
form  as 

lyiw; 


Wn(w)r, 


(Ja) 


W) 

or  Wjj.(w)= _ (3b/ 

<L  ^(w)^?(w)  +  ^(w^^^CwJ+a2  ^(w)],-  a2^(w) 

With  a  finite  number  of  samples,  one  spectral  density  estimate  is  the 
periodogram.  A  good  spectral  estimate  is  obtained  by  averaging  several 
periodogram  from  segments,  denoted  as  M,  that  provides  the  best  spectral 
density  estimate  and  the  time- delay  estimates 

for  deterministic  signals, the  resulting  optimum  filters  are  in 
terms  of  signal  and  noise  spectral  densities. 


WI(w)= _ 

in1(u)5U(") 


5y(w) 

+  ^(w)(f  ^(wJ+a^^Cw)) 


(4a) 


or 


W  (w)»_ 


£y(w; 


5^  (w)lz2(w)-a2^(w) 


(4b) 


the  difference  between  W^.  and  is  ^ue  £0  variance  of  the  signal 

which  is  negligible  for  small  a'  and/or  low  signal-to-nolse  spectra* 


Otherwise,  both  filters  are  similar  where  the  firnt  term  within  brackets 
in  the  denominator  is  due  to  ,,(w)  and  the  second  term  <j)  (w).  Hasoab  and 

Boucher[4 3  have  examined  the  behavior  of  every  window  in  the  presence  of  a  strong 
spectral  peak  in  $  (w).  With  a  no  noise  case  W^(w)  is  undefined  and  Wj.^(w)  is 

y 

The  expected  generlized  correlator  output  becomes 
Z^(w)Z»^(w)=  a^(w)e"^WC” 

divided  by  (w)( since  the  window  W  T(w0  =  "— — — )  then  ,  we  get  the  term  of  e  ^w^", 

*y  ^  (w) 

the  time  delay  estimation.  Hence  in  this  caseywe  use  W^Cw)  to  process  the  signal 
of  sinusdid  jin’  - 

In  the  case  of  signal  unknown  Burg’s  method  is  U6ed  to  estimate  the 

J  „(w)  and  <£z2(w). 

Xz1 


reduced  to  1  _ 

§Tv) 


MAXIMUM  ENTROPY  SPECTRUM  ESTIMATION:  Burg's  Algorithm 


The  maximum  entropy  mnthod(MEM)  using  the  Burg's  algorithm  was  applied  by 
Ulrych^]  ,who  showed  the  remarable  resolution  properties  of  this  approach.  In 
applying  the  concept  of  maximum  entropy  to  spectral  analysis  (6]  we  begin  with  the 
relationship  between  the  entropy (strictly  speaking,  the  entropy  rate  for  an  infinite 
proaess)  and  the  spectral  density  S(f)  of  a  stationary  Gaussian  process, 

log  S( f)df:  .  (5) 

where  f^  is  the  Nyquist  frequency. 

Rewritting  (5)  in  terms  of  the  autocorrelation  ip(k)  of  the  process  gives 


whereat  is  the  uniform  sampling  rate.  Maximizing  (6)  with  respect  to  the  unknown 
(|(k)  with  the  constraint  that  S(f)  must  also  be  consistent  with  the  known  auto¬ 
correlations  <t>(0),<t>(l results  in  the  MEM  spectral  estimate. 


This  estimate  expresses  maximum  uncertainty  with  respect  to  the  unknown 


information  but  is  consistent  with  the  known  information.  The  variational  procedure 


The  estimate  expresses  maximum  uncertainty  with  respect  to  the  unknown  information 
but  is  consistent  with  the  known  information.  The  variational  procedure  leads  to  the  r. 
expression  for  the  spectral  densi ty[_6  ],  [_?]*  which  for  n  real  linear  process  x(i) 


PE(f)-* 


l+V^r^  exp(-j2nfiAt)J  2 

7*\ 


(?) 


where  PM  is  a  constant  and  r^  are  prediction  error  coefficients  that  are  determined 
from  the  data. 


The  chief  shortcoming  of  the  MEM  spectral  estimates  has  been  the  lack  of  a 
quantitative  method  of  determing  the  length  of  the  prediction  error  filter  r(t) 
in  (7).  Recent  work  of  AR  processes  appears  to  overcome  this  problem. 


The  algorithm  of  Burg's  method  is  considered  a  wide-sense  stationary  Gaussian 
process  x(t)  of  zero  mean  and  duration  T  seconds.  The  detail  of  this  algorithm  is  in 
reference  _ 


III.  Experimental  Study 
(A)  Signal  known? 

The*. random  signal  is  generated  by  passing  white  noise  through  a  recursive 
filter  with  the  impulse  response  Bnexp(-Bn).  The  Gaussian  white  noise  at  the  outputs 
of  two  sensors  are  independent  and  have  been  generated  by  Monte  Carlo's  method.  An 
implementation  of  the  generlized  correlator  is  executed  fbr  each  set  of  128  data 
points  obtained  from  the  sensor  outputs.  The  time  delay  between  the  sensor  elements 
is  arbitrary  and'  is  selected  as  16  &  10  units  with  di  fferent  case.  Here  to  aompare 
the  time  delay  results  from  different  signal-to-noise  ratios  and  to  determine  the 
optimum  M. 

When  the  signal  y(t)  is  random  with  known  statistics  plus  the  Gaussian 
noise  with  zero  mean  and  variance  ,  the  signnl-to-noi sa  ratios  of  -10dB,0dB,and 
10dB  are  chosen  with  results  given  in  Figs.  2-6  with  B=3. 55,r=l6.  Figs. 7-11  with 

B»17,T*10. 

The  sinusoidal  signal  y(t)  is  known  and  plus  the  Gaussian  noise,  we  have 
the  results  shown  in  Figs.  1?-  |jj. 


(B)  Signal  unknown! 

From  the  parameter  estimator^]  which  with  the  addition  of  noiRe  in 
the  output  of  the  AR  model. is 

& ( n ) -y ( n )  +  n ( n ) 

Hence,  the  model  for  the  sinusoids  in  white  noise  problem  contains 
poles  and  zeros.  In  using  AR  parameter  estimation  methods  that  assume  all¬ 
pole  model  is  dif  ficnlty « ('?]  [lOj  It  is  well  known  that  the  least  squares  <- 
estimateof  coefficient  is  biased  in  this  case,  and  the  bins  is  caused  by 

.  m  . 

the  expected  value  of  1  i  0  (  W-a/i>  £  ).  Correspond  i  ng  to  the  FPE  we  can 

get  <T(w)  hy  substracting  ^(w)  with  FPE. 


'J(n) 

_ •>— 

y(n) 

White 

noise 

Wavelet 

1  1 

(a)  AR  model 


Blo^k  diagram  (a)  AR  model  (b)  ARMA  Model 

Then,  use  this  power  spectrum  as  the  window  parameter  of  Hassab  & 
Boucher  optimum  window;,  and  we  can  estimate  the  time  delay,  from  a  generlized 
correlator. 

Similar  to  the  known  signal  case  when  y(t)  in  random  with  iitiknown 
statistics  plus  the  Gaussian  noise  with  zero  mean  and  variance,  the  signal-to- 
noise  ratios  of  -10d3,OdB,and  lOiiB  are  chosen  with  results  given  in  Figs.  (5-/^ 
with  B=!?.33,'E*=l6.  In  Figs  20  lfwi  th  B-1?  ,T=10. 

The  sinusoidal  signal  y(t)  is  unknown  and  plus  the  unknown  white 
Gaussian  noise,  we  have  t lie  results  shown  in  Figs.  25- ,withc'-l6. 


IV.  Summary 


The  time  delay  estimation  from  cross-correlator  without  window  exhibits 
a  thresholding  effect  as  the  eignal-to-noise  ratio  decreases.  In  which  we  get  the 
false  peak. ^11 ^ 


After  using  the  window  and  averaging  the  periodgrams  of  the  data  segments 
the  thresholding  effect  can  be  avoided  then  the  maximum  M  is  chose.  The  unknown 
signal  case  before  we  compute  the  (£^(w)  the  Akaike  FPE  criterion  must  be  used  to 
decide  the  order  of  the  prediction  error  filter  coefficients. (Seo  Table  1.) 

The  relation  of  SNR  and  M  can  be  found  in  Fig,^$  from  the  Fig.  we  found 
M  is  increased  when  SNR  is  decreased. 

In  the  random  signal  case  and  11=17  are  chose*  Here  we  can  find 

B=1?  has  the  small  power  comnare<i  to  and  the  small  M ( =7 )  is  enough  to 

detect  the  time-delnv  with  SNR--30dB. 


The  key  of  time  delay  estimation  of  tinknown  signal  is  the  estimation  of 
$  (w),  and  $  ^(w).  If  we  can  estimate  the  power  spectrum  of  y(t)  from 

the  sensor:  outputs ,  than  we  must  have  a  good  estimate  using  the  above  method. 
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1  Block  diagram  of  the  generalized  cross-correlator. 

?n  Without  window  output  for  M-1,  -16,  3NR=10dB, the  result  of  time-delay  =16. 

2b  window  function  of  W  (w),  for  M=1,'t=t6,  SNR=lodb. 

2c  With  window  of  Wj - (wJ ,M=1 ,BNR=10dB,  =16,  the  result  of  time-delay  =16. 

3a  Without  window  output  forT=l6,  SNR=OdB,  M=1 ,  the  result  of  time-delay  =10. 

3b  Window  function  of  W  _(w),  for  M=1,T=16,  3NR=0dB. 

3c  With  window  of  W^^wy,  for  M=1,T=l6,  SNR=OdB,the  result  of  time-delay  =13. 
4e  The  same  as  Fig. 3c, but  M=3,  the  result  of  time-delay -» 13. 

4b  The  same  as  Fig.  ^c,  but  M=7,  the  result  of  time-delay  =16. 

5a  Without  window  output  for  =16,M=1 ,SNR=_iOdB,  the  result  of  time-delay  =-6. 
3b  Window  function  of  W  (w),  for  M=1,*t=l6,  SNR=-10dB. 

5c  With  window  of  W  (v),  M-1  ,SNR=-10dB,  1  =  16,  the  result  ofT=15. 

6a  With  :  window  out  out  for  1-1=7,  BNR=-lOdB,  the  result  of  time-delay  =12. 

6b  With  window  for  M=15,  SNR=-10dB,  the  result  of  time-delay  =13. 

6c  With  window  for  M=20,  SNR=-10dB,  the  result  of  time-delay  *l4, 

7a  Without  window  B=17.  T  =  10,  .SNR=10dB,M=1 ,  the  result  of  time-delay  =10. 

7b  With  window  B-17.,  1  =1o1SNR-10,M=1 , the  result  of  time  delay  =10. 

8a  Without  window,  for  B.=  17,  oNR~0rtB,  =10,M=1,the  result  of  time-delay  =19. 

8b  •  <  i  With  window  for  B- 1 ? . ,  BNR-OdB, T  =10,M=1 , the  result  of  time  delay  =10. 

9a  Without  window  , for  B=17,  SNK=-10dB,  =10,  the  result  of  time-delay  =13* 

9b  With  window  for  B--17,  SNR=-1()dB,  =10,  M=1,  the  result  of  time-delay  =10, 

10a  Without  window  B=1V>  0NRr--2C>db,M=1 ,  =10,  the  result  of  time-delay  =13* 

10b  With  window  B=17,  8NR=-20dB,  M=i . T  =  10,the  result  of  time-delay  =13. 

10c  With  window  B=17,  SNR=-20dB,  M=5, T=1°.  the  result  of  time  delay  =10. 

11a  Without  window  B=17,  SNH=-30dB,  M=7,  *t  =  10,  the  result  of  time  delay  =22. 

12a  Without  window  sine  wave  for  M=1 ,SNR=5dB,  ?  =  l6,  the  result  of  time  delay  =16. 
d2b  With  window  sine  wave  for  M=1,  KNR=9dB, T =16,  the  result  of  time  delay  =16. 
13a1  Without  window  sine  wave  for  M=1,  SNR=-5dD,T=l6,  the  result  of  time  delay  =11 
13b  With  window  of  sine  wave  for  M=1 ,  SNR=-'jdB,?=l6,  the  result  of  time  delay  =1? 
14a  Without  window  of  sine  wave  for  M=1 ,SNR=-15dB, T =16,  the  result  of  time  delay 
=-20. 

14o  With  window  of  sine  wave  for  M=20,  BNR=-15dB, "f=l6,  the  result  of  time  delay 

=  12. 

19a  Without  window  output  signal  unknown  for  B=3*33,M=1 ,7=l6,BNR=10d&,  the  time- 
delay  =  16. 

15b  Window  function  of  W .,(w),  for  B=3«33.8HH=10db, T*l6« 

19c  With  window  of  Wj._.(wV,  for  B=3.33,SNR=10dB,  T‘s'6,  the  time  delay  =16. 
l6a  Without  window  output -(if  B=3.33,  SNR=0dB,  T  =  1 8, (unknown  signal),  the  result 
of  time  delay  =13.  M=1. 

16b  Window  function  of  W  (w),for  B=3,35,  SNR-;0dB,  J=l6,  (unknown  signal) ,' Mai . 

16c  With  window  of  W  (w;  for  B=3.33»  3NR=0dB,  T~16,  (unknown  signal),  the  result 
of  time  delay  =  13  ,M-1. 

17a  Without  window  output  with  unknown  signal,  for  B=3. 33,SNR=OdB,  T=l6,  the 
result  of  time  delay  =15,M=3. 

17b  The  same  as  17. a  but  M=H,  the  result  of  time  delay  =15. 


Fig. 

Fig. 


Fig. 

Fig. 


17c  The  same  as  17. abb, but  M=7,  the  result  of  time  delay  =l4. 

l8a  Without  window  of  unknown  signal  Bt3.33»SNR=-10dB, rt=l6,  M=1,  the  result  of 
time  delay  =-6. 

l8b  Window  function  of  W  (w)  , B=3«93,  CNH=-10dB,  =l6,M=1,  unknown  signal. 
l8c  With  window  of  W^(w7,  B=3. 33,8NR=-10dB,  =  l6,M*T,  the  result  of  time  delay 


=11. 
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Fig. 

26b 

Fig. 
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Fig. 

27b 

Fig. 

28 

Note: 

Fig.O 

With  window  ,the  same  as  18-c,  but  M=7 ,  the  result  of  time  delay ^=11 • 

Tho  seme  ns  19-a,  but  H=15,  the  result  of  time  delay?=l4. 

Without  window  of  unknown  signal  for  Bt17.  ,3NR=10dB,  M=1 ,  the  result 

of  time  delay7  =  10. 

With  window  of  unknown  signal  for  B=17,SNR=10dB,‘Ts10,M=1 ,  the  result  of 
time  delay ^£=10. 

Without  window  of  unknown  signal,  for  13=17,  SNR=odB,^=lO,  KM,  the  result 
of  time  delay  <T  =  10. 

With  window  of  unknown  signal,  for  13=17,  SNR=Odb,  "J  =  10,M=1 ,  the  result  of 
time  delay  T=10. 

Without  window  of  unknown  signal,  for  B=1?.,  SNR=-10dB,  M=1,  the  result 
of  time  delayT-13.  — 

With  window  of  unknown  signal,  for  B=17,  SNR=-10dB,  (  =10. ,  the  result  of 
time  delay 7~10,  M=3. 


Without  window  of  unknown  signal,  for  B=17,  SNR=-?0dB,  "^=10,M=1 ,  the  result 
of  time  delay  ^  =  15. 

With  window  of  unknown  signal,  for  B=17,  UNR=-20dB,  =10,M=1,  the  result  of 
time  delay  =10. 


Without  window  of  unknown  signal  for  B=17,  SNR=-30dB, 7”=10,  Hal,  the  result 
of  time  delay 7  =-6. 

Without  window  of  unknown  signal  for  B=17,  SNR=-30dB,  TMO,  Ma?,  the  result 
of  time  delay  7=22. 

With  window  of  unknown  signal  for  B=17,  SNR=-30dB,  ^=10,  M=7,  the  result  of 
time  delnyT=1?. 


Without  window  of  sine  wave  of  unknown  signal,  for"7=l6,M=1 ,  the  result  of 
time  delav~£=l6.  SNR=cidB. 

With  window  of  sine  wave  of  unknown  signal,  forV=l6,  Hal,  the  result  of 
time  delay  =  16.  3NR=9dB, 


With  window  of  unknown  signal  M=1,  the  result  of  time  delayf  =15. 
(without  window  case  is  the  same  as  Fig, 17a.)  (3NR=-SdB). 

With  window  M=3,  the  result  of  time  delay falfi.  (SNR=-3dB) 

With  window  for  SNH=-15db,7  =  1(>,M=3,  the  result  of  time  delayT=l4. 

Without  window  for  SNR=-19dB,M=3,  *16,  the  result  of  time  delay  T=12. 

The  relationship  of  random  signal  case(btexn(-Bt ). (known  signal) 


is  the  random  signal  case' with  B=3. 33  and  signal  is  known. 


'T3 


Reference  A 


(A)  Burg* a  Recursive  Algorithm 


Let  X(t)  denote  a  sample  function  of  this  process,  which  is  sampled 
uniformly  at  the  rate  l/T  to  produce  a  time  series  consisting  of  N  samples, as 
shown  by 

£x(n)J=  jx(l),X(2), . ,X(N  )J  (Al) 

where  X(n)  denotes  the  nth  sample  of  the  series. 


The  autocorrelation  function  of  the  time  series  X(n)  for  a  time 

lag  of  mT  seconds  and  time  index  n  is  defined  by 
s 

Rx(n,m)=  K  (X(n,m)  X*(n)j  (A?) 

where  E(_  ‘Jdenotes  the  expectation,  and  the  asterisk  denotes  the  complex 
conjugate  operation,  With  the  random  process  X(t)  assumed  wide-sonse  stationary. 


the  autocorrelation  function  R  (n,m)  dependent  only  on  the  time  lag  mT  ,  as  shown 

A  R 

by 


R^(n,m)=H^(m)  (A3) 

Also,  the  autocorrelation  function  exhibits  conjugate  symmetry,  i.,  R^(m)  is 
a  Hermitian  function 


* 

Rv(-m)=K*(m)  (Ah) 

Suppose  that  the  time  series  |x(n)^  is  applied  to  a  linear  digital  filter  of 
impulse  response  ^h(n)j  of  order  P.  The  filter  produces  an  output  time  series 
|  ,y.(n)  ^  that  is  designed  to  approximate  a  desired  time  series  |  d(n)|  ns  shown 
in  Figure  Al.  and  the  resultant  error  time  series  is 

e(n)  =  d(n)-.y(n)  ( A*? ) 


X(n) 


Y(n) 


d(n) 


e(n) 


Fig.  Al 


-7 


For  a  special  case  of  a  filter  designed  to  predict  the  value  of  a 
random  process  X(t),  one  time  unit  ahead,  by  using  the  present  and  past  values 
of  the  time  series  X(n),  i.e. 

d(n)=X(n+1 ) 

A  prediction-error  filter  is  introduced.  Fig.  A2  shows  the  functional 
relationship  between  the  linear  prediction  filter,  characterized  by  the  impulse 


Fig.A2  :  Illustrating  the  relationship  between 

the  predictive  filter  and  prediction-error  filter 


The  prediction-error  filter  aquation  expressed  in  matrix  form  is  given 

by 


(A6) 


where 


M-i 

r  12 
K-o 


(K) 


V-K) 


is  the  prediction-error  power 


A  recursive  algorithm  developed  by  Burg  ^ I}"]  for  solving  the  set  of  prediction- 
error  filter  equation  io  availible. 


Suppose  that  we  know  the  solution  to  the  set  of  p  equations  pertaining  to 
a  prediction  filter  of  order  p,  as  shown  by  (46).  We  can  take  the  complex  conjugate 
of  both  sides  of  (A6)  and  recognize  that 


(1)  The  prediction-error  power  is  a  real  quantity 

(2)  The  autocorrelation  matrix  is  Hermitian,i.e. 


R^  *  ( m )  =  R^  ( -m ) 

We  thus  obtain 


("  V°>  V1}  . 

Hx < “ 1 }  Kx(0)  . 

...  R  (M-l) 
...  RX(H-P) 

Vi(1) 

• 

• 

— 

'  l  > 

*M 

0 

* 

• 

« 

H,  (1-M)RV(2-M).:... 
t*  a  X 

...  Hx(0) 

1 

w..  (  M- 1 ) 

v  ,l1-1 

♦ 

.  0  . 

(A  7)  can  be  rearranged  as 


^  V 

(  Rx(0)  Hx(-1)  .  RX(1-H) 

Wp_  ^ ( M— 1 ) 

0 

fix(l)  Rx(0)  .  Rx(2-M) 

M-1 

0 

♦ 

• 

— 

9 

r 

I 

RX(M-1)  KX(K-?) .  Rx(0) 

1 

U ) 

(A6)  pertains  to  a  prediction-error  filter  of  order  p  operated  (P~  m) 


(A7) 


(Aft) 


in  the  forward  direction.  On  the  other  hand,  (a8)  pertainn  to  a  prediction- 
error  filter  of  the  name  order  except  that  it?  operated  in  the  backward 
direction.  This  is  illustrated  in  t'ig.  A3 


|  X(1),X(^) 


i'ig.A3:  illustrating  ;  (a)  forward  and  (h)  hackwarrt 
prediction  error  filtering 

Combining  (Ab)  and  (AS)  to  expand  the  number  of  prediction  error  filter 
equations  by  one  is  as  follows: 


r Hy(0) 

V-1' . 

V°>  . 

K  ( 1  -!•: ) 

aA?-k> 

Rx(-M) 

R..(1-M) 

A 

s 

s  'b 

1 

W  ( •’ 

M_  1 

• 

•  ' 
• 

+  IV 

'  N,\ 

0 

W  \.  (l) 

r  -1 

•  1 

^RX(K.) 

H  (h-,P) . 

RX(M-D . 

ky(o) 

s«> 

Rx(-1)  j 
Rx(°)  J 

w.,  (M- 

h-1 

0 

\ 

1 ) 

( p  ) 

r  *  "i 

N 

M 

h 

0 

0 

=  1 

+  Wm(M) 

However,  for  the  correspond inp  prediction-error  filter  of  order  (m+l)  we 
have 


1  R x(0) 

R^-l)  .  .. 

P,{(1-H) 

V 

Rx(-M) 

'  1  ' 

'  \ 
P 

M+1 

HX(l) 

Rx(0)  . 

Kx(2-N) 

RX(1-M 

V" 

• 

0 

• 

R.(M-I) 

X 

RX(M-P) _ 

Rx(0) 

• 

Wn(M-1) 

• 

• 

0 

RX(M) 

. 

HX(D 

Ry(0) 

I 

w  (M) 

k  J 

0 

(AIO) 


Uompnrinp  with  (Ay)  and  (AIO),  we  deduce  that 


W„(K)=  W  (K)  +  W.  (M)  W.,  (M-K) 

n  M-1  M  M-1 

;  K=0, 1 , . . . . ,  p(-sli ) 

(All) 

P  -  P  4  W  (M)A  * 

M+1  M  M  M 

(A12) 

and 

0=-4*  +  W  ( h)  P,, 

M  h  M 

(A13) 

In  the  recursive  formula  of  (All),  note  that 

for  all  values  of  M 

(  1  for  K=0 

VK)  ~  i 

(  0  for  K>M 

(A13)  will  always  has  a  solution  provided  that  P^  >0,  thus  us.inp  (Al$)  to 
eliminate  A*  from  (Al?),we  pet 


By  analopy  with  the  transmission  of  power  through  a  terminated  two-port 
network,  we  may  view  W  (M)hb  a  "  reflection  coefficient". 


COMPUTATIONAL  PHUCi'.DU H b 


Summarizing  the  steps  in  the  recursive  procedure  involved  in 
calculating  the  prediction-error  filter  coefficients  and  related  values 
of  the  spectral  density  and  autocorrelation  function,  they  are 

(l)  For  thegiven  time  series  X(n)  ,  n=1,P, . N  and  Mr-O. 

Calculate  the  power 

N 

P1=T*Zx(n)  X* ( n) 
n=\ 

and  calculate  the  reflection  coefficient  and  the  prediction-error 


power  by 


W  = 


X*  (n)  X(n+l) 

*  -  \ 


L[|X(n)|  +  |x(n+l)|  ^ 

(?)  Increament  M  by  1,  and  calculate  the  next  values: 

Vn)  =  PM-1(n)  +  WH(M)  qM-1(n) 

qM(n)=qM-1(nM)  +  WM(M)  ?M-1(n+l) 


where 


P^Cn)*  X(n) 
q^(n)=  X(n+l) 


(3)  Calculate  the  refle  ju  coefficient 


N-  M 

XjPM-l(n) 


W  (M)= _ — 

V’  N-M 


22  j|PM-1(n)l'J  *  |qM-1(n)|  ?  1 

-i/i  -  t  * 


and  the  corresponding  value  of  the  error  prediction  power 


W  h,  MVH)I2 1 


(4)  Calculate  the  remaining  coefficients  of  the  prediction  error  filter 
by 

W„(K)  =  WfU1(K)  +  Wm(M)  WH-i(M-K)  K=1,? . K-1 

(5)  Repo&t  steps  (?)  to  (4)  for  each  value  of  M  up  to  the  optimum  value 
p(=M)  for  which  the  final  aarediction-error 


(l-UM+i ) 


( N-H-1 ) 


M+1 


as  a  minimum. 


(6)  We  now  have  all  the  quantities  required  for  calcuation  tne  extrapolated 
value  of  the  autocorrelation  function. 

KM 

KX(K)  -  PM  -  £  WH  (K)  HX(K-K) 

‘^  =  1 

and  for  calculating  the  spectral  density  estimate 

P. 


M+1 


2B  1  1+ ^  WM(K)  exp(-?lt  kifp)  2 

*  fc*  \ 


where  T  is  the  sampling  period  and  X(t)  is  the  band-limited  to  -B,B. 
6 


When  the  process  X(t)  is  sampled  at  the  Nyquist  rate,  we  have  T=l/?B 
The  computation  is  thereby  completed. 
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